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ABSTRACT 

By using various properties of the complete elliptic integrals, we have derived 
an alternative expression for the gravitational potential of axially symmetric bodies, 
which is free of singular kernel in contrast with the classical form. This is mainly a 
radial integral of the local surface density weighted by a regular "mean Green func- 
tion" which depends explicitly on the body's vertical thickness. Rigorously, this result 
stands for a wide variety of configurations, as soon as the density structure is vertically 
homogeneous. Nevertheless, the sensitivity to vertical stratification — the Gaussian 
profile has been considered — appears weak provided that the surface density is con- 
served. For bodies with small aspect ratio (i.e. geometrically thin discs), a first-order 
Taylor expansion furnishes an excellent approximation for this mean Green function, 
the absolute error being of the fourth order in the aspect ratio. This formula is there- 
fore well suited to studying the structure of self-gravitating discs and rings in the spirit 
of the "standard model of thin discs" where the vertical structure is often ignored ,but 
it remains accurate for discs and tori of finite thickness. This approximation which 
perfectly saves the properties of Newton's law everywhere (in particular at large sep- 
arations), is also very useful for dynamical studies where the body is just a source of 
gravity acting on external test particles. 
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1 INTRODUCTION 

The capability to calculate properly the gravitational po- 
tential inside celestial bodies is a longstanding challenge 
in Astrophysics, mainly due to the difficulty to manage 
the hyperbolic divergence l/\r — f\ (a major feature of 
Newton's law). The presence of strong inhomogeneities of 
density and complex geometries or shapes seem less prob- 
lematic in comparison. For a few special configurations 
(not always realis tic), the potential/density pair is sim - 
ple and analytical (|Mestellll963l : iBinnev fc Tremaine!ll987l ). 
However, for most realistic matter distributions, a nu- 
merical computation of the gravitational potential has to 
be performed, ei ther from the integral form when accu- 
racy is required ( Ambastha fc Varmalll985 : iHachisul 1 19861 : 



ISto'ne fc NormarJll992t iBaruteau fc Massed |200S|) . or from 
the Poisso n equation, when computational time is th e limit- 
ing factor |Bodo fc Curirll 19921 : [Fromang et al.ll2004) . While 
most theoretical developments allow to reach any degree of 
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accuracy, their numerical implementation are usually limited 
by errors and uncertainties of various origins. A striking il- 
lustration is the infinite series representation of the Green 
function in terms of Legendre polynomials which, although 
exact, is known to h ave a (very) low convergence rate in- 
side sources (see e.g. iKellogd Il929l : | Durandl 1 19531): so. its 
implementation is o ften disappointing in practice (|Clementl 
11974 lHachisulll986D . Compact integral expressions, despite 
singularities, may the n be preferred l|Cohl fc Tohline|[l999l : 
iHure fc Pierensll2005D . 

The gravitational potential of discs, starting with galac- 
tic discs, has received much attention in the last decades 
(|Mestelll 19631 : IBinnev fc Tremaindll987l ) . Even for very sim- 
ple and academic configurations, new analyti c solutions and 
appro ximat ions co ntinue to be pr oduce d (ICox fc Gomez 



Hure et al 



| 2002 j: iSchulzl 120091: IVogt fc Letelier|l20iq , 12011 
hoill; iBannikova et al II2OIII : ISchuhjhoill ). either for estab- 
lishing theoretical references, or for analysing observations. 
The problem is far from being solved. Searching for new 
formulae or techniques aiming at improving the description 
(in terms of computation time, accuracy, series convergence, 
singularity avoidance, etc.) of gravitational potentials and 
forces of celestial bodies is therefore still of relevance today. 
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Figure 1. Typical configuration for the gravitating, axially sym- ' mo dulus k 

metric body (finite size and mass, and local total thickness 2h) 

and notations associated with the cylindrical coordinate system. Figure 2. Variation of feK(fc) and fcE(fe) with the modulus k. The 

first function logarithmically diverges as k — > 1, while the second 
function is bounded (and takes its maximal value for k S3 0.91). 
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Hure fc Hersanj 12011 ). we have mainly investigated 
fiat, axially symmetric discs where the surface density 
is a power-law of the equatorial radius, with a special 
att ention to edge e ffects . In this paper, we consider, as 
in IPierens fc Hurel (|2005h ■ systems with finite thickness, 
and present a new decomposition of the gravitational 
Green kernel which permits to express the potential as 
a sum of two regular integrals. These two integrals have 
interesting properties and their expression can give rise to 
very accurate approximations in the limit of geometrically 
thin discs, i.e. when the vertical thickness h of the disc is 
small enough. 

This paper is organized as follows. In Section [2] we re- 
call the general expression for the potential valid under axial 
symmetry and define the associated "mean Green function" 
to be calculated. In Section [3] we demonstrate that this 
mean function — an improper integral — can be rewrit- 
ten as the sum of a regular, line integral over the boundary 
of the system and of an analytical function, also regular. 
This enables to rewrite the expression for the potential in a 
different, more tractable form which contains a surface in- 
tegral and a line integral. This is the aim of Section [4] In 
Section [5] we briefly verify that the new expression for the 
potential automatically fulfills the right conditions at infin- 
ity. The question of the reduction of the remaining surface 
integral is addressed in Section [6] In Section [7] we build 
an approximation for the mean Green function from a first- 
order Taylor expansion (systems with small aspect ratios). 
We check the accuracy of the approximation and associated 
potential in Section [8] through a typical toroidal configura- 
tion. In Section[9j we show that stratification effects remain 
weak inside and outside the source, as soon as the total sur- 
face density is conserved. The last Section is devoted to a 
conclusion. A few useful formulae and demonstrations are 
found in Appendix. 



2 THE GRAVITATIONAL POTENTIAL 
UNDER AXIAL SYMMETRY. 
THEORETICAL GROUNDS 

We consider an axially symmetric disc defined by its merid- 
ional cross-section (fi) and mass density p, as depicted in 
Fig. [1] The gravitational potential ip of this body is given, 
in cylindrical coordinates (R, Z), by the surface integral (e.g. 
lDurand|[l95l ): 



i>(R,Z) 



zp(a, z)kK(k)dadz, (1) 



where (a, z) refers to points inside (fi) belonging to the 
source, K(fe) is the complete elliptical integral of the first 
kind (see the Appendix [X] for its definition), 



2VaR 



y/(a + R)* + C- 



€ [0, 1], 



(2) 



is the modulus, and £ = Z — z is the altitude difference. 
The divergence of the term l/\r — f*| when r — > r initially 
present in Newton's law still exists here through the func- 
tion K(fc) which logarithmically diverge^] as k — > 1 (corre- 
sponding to a = R and z — Z). This singularity is clearly 
visible in Fig. [2] where we have plotted the function feK(fc) 
versus k. Although K(fc) — >■ oo as k — > 1, the i ntegral of 
the potential is finite (|Kellogdl 19291 ; |Purandlll953l ). Outside 
the material domain (fi), Eq.fTJ) can be easily evaluated nu- 
merically as the modulus k is always less than unity. But 
difficulties begin as one approaches the boundary (dQ). In- 
side the source, the direct computation of the integral by 
standard numerical techniques cannot lead to accurate po- 
tential values due to the improper integral (at the location 
where k reaches unity). The separate treatment of the log- 
arithmic divergence (see note[T]) gives good results in terms 



Precisely, we have l lGradshtevn fc Rvzhikll 19651) : 
K(fc) ~ In4-In>/1- fc 2 



(3) 
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of pre cision (|Ansorg et al . 2003; Hure 2005; Bannikova et al.l 
l201ll ). but this approach is not very tractable. For these rea- 
sons, the Green function is generally replaced by its series 
expansion in Legendre polynomials. This effectively removes 
the divergence but generates new problems in practice (low 
converg ence rate of the alternate se ries, truncations, etc.; 
see e.g. IClemenj[l97i lHachisdli986h . 

The way to perform the double integral in Eq.fTJ de- 
pends more or less on the mass density distribution and 
body's shape via the equation for its boundary (dQ). If this 
equation can be defined by two bijections of the form z± (a) 
(i.e. z+ above the equatorial plane and Z- below), then we 
can perform the integral over z first, followed by the integral 
over a. This is the most natural approach to treat geomet- 
rically thin d iscs and rings (e.g. IShakura fc Sunvaevlll973l ; 
|Prinridfl98ll ). and this is the case we will consider in the 
following. The body is then regarded as a collection of in- 
finitely thin cylinders. Assuming that the mass density does 
not depend upon z but only varies with the radius a, Eq.{l]) 
can be written in the form: 



ip(R,Z) = -G I g (R,Z;a,z±(a))Z(a)da, 

J a 

where 

E(a) = / p(a)dz 
J 2 

is the total surface density in the disc, and 



«—kK{k)dz 



dz 



(4) 



(5) 



(6) 



plays the role of a "mean Green function". Actually, 
K(fe) rises around a = R and £ = 0, while it has small 
amplitude elsewhere (K(fc) — > ^ for small fc). Due to the 
integration process, Q is expected to be a regular function, 
peaking around a = R (see below). It seems that this integral 
is not known in closed-form, but only through an alternate 
series whose convergence rate is unfortunately low IjDurandl 
[195M 

For certain configurations where z±(a) is not bijective, 
it may be more convenient to integrate over a first, and then 
over z. This procedure is equivalent to considering the body 
as infinitely flat discs piled up along the z-directiorjf]. Fi- 
nally, it is worth mentioning that the double integral over 
the meridional cross-section can be converted, under cer- 
tain conditions, into a line integral over the boundary (cT2) 
through the curl theorem which writes: 



(d a N -d z M)dadz= / (Mda + Ndz), (7) 

'(f2) J(dCl) 

under axial symmetry. Here, this requires the existence of 
two cylindrical functions M(a, z) and N(a,z) such that: 



- 2, 



-pkK(k) = d a N - d z M. 



(8) 



2 Actually, in contrast with the two components of the gravita- 
tional acceleration due to a thin cylind er, the potenti al is appar- 
ently available only through the series l lDuran d 1953). 

3 The f o rmula for homo geneous, infinitely th i n disc is given in 
Durandl dl953ll (see also lLass fc Blitzed Il983l : iFukushimal l20ld : 
Tresaco et alj|201lfl . 



This is not guaranteed in the general case where d ensity 
gradients are present. However, lAnsorg et al.l ((2003) have 
solved this question in the fully homogeneous case, i.e. when 
p is a constant. 



3 BYPASSING THE KERNEL SINGULARITY: 
THE CASE OF VERTICALLY 
HOMOGENEOUS SYSTEMS 

The absence of any closed-form for Q is problematic for 
most theoretical and numerical applications. We can how- 
ever rewrite this mean Green kernel in another, more conve- 
nient form, as follows. By using various relationships involv- 
ing the complete elliptic integrals of the first, second and 
third kind (K, E and n respectively; see the Appendix |A"1 
for their definition), we have established by direct calculus 
the following two formulae, valid for any integers n and p: 

f- n - 1 uP 

d c Ck p K(k) = i — {[(n-p+ l)m 2 + (p - l)fc 2 ] K(k) 



, 2 2 , E(fc) 

+ (fc - m — V 1 



(9) 



and 



d ( Ck p n{m,k) = 



+k 



m 
2 E(fc) 



{ [(n - p)m 2 + (p- l)fc 2 ] n(m, fc) 

(10) 



where 



2VaR 



a + R' 

is the characteristic or parameter of n, and 



fc' = s/l - fc 2 



(11) 



(12) 



is the modulus complementary of fc. Note that m — > fc as 
C, — > and ^ ^ m ^ 1. If we now combine together 
Eqs.© and (llOf) . with weights 1 and m' 2 = 1 — m 2 respec- 
tively, we get the general formula: 



a c C n fc p K(fc) - m' 2 d c Ck p n(m,k) 



C n_1 fc p 



(13) 



x { [(n - p + l)m 2 + (p - l)fc 2 ] K(fc) 

-m' 2 [(n - p)m 2 + (p - l)fc 2 ] n(m, fc) - m 2 E(fc) j . 

The remarkable point is that, for p = 1 and n — p, the 
complete elliptic of the third kind disappears in the right- 
hand-side, and we get: 

<9 c CfcK(fc) - m' 2 <9 c <fcn(fc) = fcK(fc) - fcE(fc). (14) 

As m does not dependent on the relative altitude we can 
write this relation as: 

c\ {< [fcK(fc) - m' 2 fcn(fc)] j = fcK(fc) - fcE(fc). (15) 
Finally, given d£ = — dz, it follows tha10: 

kK(k)dz= [ kE{k)dz-(U(m,k) (16) 



4 iHure fc Dieckmannl d2012f ) make another use of this relation- 
ship. 
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where H is defined by 

H(m, k) — k [K(fe) - m' 2 n(m, fc)] . (17) 



the well-known expression: 



We recognize in the left-hand-side of Eq. (fl6|) the ex- 
pression coming into the definition of our "mean Green func- 
tion", in Eq.©. This new form is very interesting. Actually, 
the first term in the right-hand-side of Eq. (|16[) is fully reg- 
ular since E(fe) £ [f,l]- The function fcE(fc) is displayed 
versus k in Fig. and it is to be compared with the ini- 
tial kernel fcK(fc). The second term is also fully regular and 
bounded. The divergence of H(m,fc) when k and m both 
approach unity is cancelled out by the presence of the two 
vanishing factors m' 2 and f. The H-function is plotted ver- 
sus m and k/m in Appendix [B] We conclude that Eq. ()16[) 
has a three satisfactory properties: i) it no longer contains 
any singular integrand, ii) it ensures that the effect of initial 
singularity is automatically accounted for, and iii) it fully 
saves the Newtonian properties of the potential and associ- 
ated forces. 



4 A GENERAL EXPRESSION FOR THE 
NEWTONIAN POTENTIAL OF 
VERTICALLY HOMOGENEOUS BODIES 

By inserting Eq. (|16[) in Eq.©, we find that the mean Green 
function associated with vertically homogeneous systems 
writes: 



g(R,Z;a,z ± ) = M± 



where, from Eq.©: 



fc± 



kE(k)dz 
C+H(m,fc + ) + C-H(m,fc_)], 

2y/aR 



(18) 



(19) 



and £± = Z — z± is eventually a function of a (sign + 
being associated with the top boundary, and sign — is for 
the bottom). Since k± — 1 only when £± = 0, the £H terms 
in Eq. (|18[) are always finite. It means that the integration 
of the mean Green function over the radius a — see Eq.Q 
— can be performed numerically without much difficulty as 
soon as z±(a) is known. In this integration process, some 
care must be taken though. Actually, when regarded as a 
function of a, H(m, k) is not derivable at ii = a. The jump 
in the derivative can, however, be determined (see Sect. [8] 
and Appendix [Cj . Depending on (R, Z), this jump is: 

• 0, outside the body, 

• 7r, on the boundary (dil), 

• 2-7T, inside the domain (SI). 

As a result, Q is peaked at a — R, and the maximum value 

S?max. is. 



g(R,Z;R,z±) 



(20) 



h{R) 



kE(k)dz - ( + k + K(k+) + C-k-K(k-) 



g(R,Z;a,0) = 2 J-kK(k). 



(21) 



The general expression for the gravitational potential is 
then found from Eqs. (j4| and (|18p . namely: 



,-Ui.Z) = -G j E(a)i7|da 



kE(k)dz 



(22) 



+ G 



S(a)J|i [C+H(m,fc 4 



£_H(m, k-)] da. 



We recall that this expression is exact and valid for any 
point (R, Z) of space, for any surface density profile E(o) 
(no vertical density gradients), and semi-thickness h(a). It 
contains a surface integral over (fi) and a line integral over 
(dfl) through z±(a). 

If the body is symmetric with respect to the equatorial 
plane, we have z+ = — z- = h. Unfortunately, Eo. (|18|) does 
not simplify more since fc+ and k- remain different, except 
for Z = 0. This case corresponds to the potential at the 
midplane. With = — £_ = —h, the midplane mean Green 
function writes: 



g(R,0;a,h) = 2, 



and 



S?max. — 2 



h(R) Jo 



kE(k)dz + H(m,k±) 



kE{k)dz + k±K{k ± ) 



(23) 



(24) 



5 LONG-RANGE PROPERTIES 

We can check that the long-range properties of the potential 
defined by Eq.((22| are correct. Actually, far enough from 
the system, the potential is expected to tend towards the 
potential due to a central condensation, or 

lim rip = -GM, (25) 

r— ¥QG 

where r = \/ R 2 + Z 2 is the spherical radius, and 



M = 2tt / aE(a)da 



Note that, for a flat body, z± = and so one recovers 



(26) 



is the total mass. We see that k w 2V aR/r — > as r — > oo, 
and then E(fc) — > tv/2. It means that the long-range be- 
havior is ensured by the first term in the right-hand-side of 
Eq. (|22p containing the complete elliptic integral of the sec- 
ond kind E, whereas the contribution from the H-function 
is a higher-order correction (see Appendix [D] for a more de- 
tailed analysis). 



6 FULL REDUCTION TO A 

ONE-DIMENSIONAL INTEGRAL ? 

The important question we have tried to clarify concerns the 
possibility to convert the remaining double integral in the 
right-hand-side of Eq. (|22p into a line integral (by performing 
the integral over zor(). Actually, the existence of a compact 
expression for 
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/ 



kE(k)d( 



(27) 



would be helpful, since the potential ip of a tri-dimensional, 
vertically homogeneous and axially symmetric system would 
be given by a one-dimensional integral. This is the reason 
why we have explored different paths, but unsuccessfully. 
We have for instance re-written fcE(fe) in the form of an 
infinite series over fc, or re-considered E(fc) as an integral 
over <j>, followed by t erm-b y-term integration as done in 
ICviiovic fc Klinowskil (|l994l ). However, none of these two 
approaches leads, apparently, to a closed-form in the gen- 
eral case. Equation (|27[1 can also be converted into various 
equivalent forms, like: 



fcE(fc)dC = ± 



4aR 
a + R 



E(fc)dfc 



fc 2 



(28) 



where ± = Cl/C- This integral over k is known only for 
m — 1. In particular, the indefinite form has been derived by 
Dieckmann (2011; private communicatiorjfl) , but it involves 
hypergeometric series which are not easy to manipulate. The 
definite f orrrj^l is known only when the integral bounds are 
and 1 (ICviiovic fc Klinowskil 1 19941 ; ICviiovic fc Klinowskil 
Il999l : IPrudnikov et al.lll988n which corresponds to £ in the 
range ] ± oo, 0]. This is a very special situation: a and R are 
obviously different in general , and so the c a se m = 1 is too 
restrictive. Finally, following lAnsorg et al.l (|2003h . we have 
also searched for a two-component vector (M',N') whose 
curl would yield the function y^kE(k) (see Sect. [2}, but we 
failed to find a simple and convincing answer. This question 
remains therefore open. 



7 APPLICATION TO GEOMETRICALLY 

THIN DISCS AND RINGS: A USEFUL AND 
ACCURATE APPROXIMATION 

The main assumption made up to now is the absence of ver- 
tical stratification for the mass density. It is true that, except 

5 See http: //pi .physik.uni-bonn. de/~dieckman/ 

6 The definite integral can take various equivalent form, like: 



r 



kE{k)d<; : 



r u + 

JaR I E (in sech(w)) du 

J oc 

, — r v+ E(mv)dv 

/aR / ; 

Jo w' 
\faR f 1 E(tk+)dt 
k+ Jo 



(29) 



where 



and 



tyi. 

f = (a + R) sinh u, 
v = sech(«), 

t - k 



(30) 
(31) 

(32) 



which can unfortunately not be reduced more, because formulae 
linking complete elliptic integrals with different moduli, for in- 
stance E(ax) and E(x) for any real a £ [0, 1], are missing (Landen 
and Gauss transformations are useless here). 



for very academic (i.e. unrealistic) configurations, this may 
a priori not apply to astrophysical problems — however, as 
we show in Section [9j the sensitivity to stratification remains 
weak. Possible domains of interest concern cert ain poly- 
tropic fluids (e.g. lHachisulll98rj| : Ansorg et al.ll2003r ), and ge- 
ometrically thin discs and rings (|Pringlell98ll ). For thin discs 
actually, r adial gradien ts are often neglected over vertical 
gradients (|Pringlelll98lh . at least far from edges. In this con- 
text, a vertically averaged structure wi th a uniform density 
along the z-axis is commonly assumed (IShakura fc Sunvaevl 
ll973l : ICollin-Souffrin fc Dumontll990l : lDubrullell992l ). Their 
total thickness 2h(a) = z+ — z- is everywhere small com- 
pared to the radius a, which enables many developments. 
Regarding the present problem, it means that k+ and fc_ 
are always close to each other. We can, in this case, produce 
an approximation for the mean Green function (and thus 
for the potential) by estimating Eq. (|27|l , There are various 
ways to address the problem. In order to get an approxima- 
tion valid in the whole physical space (and not only in the 
disc vicinity as often considered), it is appropriate to work 
with the mean modulus k defined by: 



k = 



k+ + k- 



(33) 



In particular, k never reaches unitjQ, preventing any 
eventual divergence of the elliptic integrals. If we now per- 
form a Taylor expansion of fcE(fc) around k, we obtain, at 
first order: 



fcE(fe) w fcE(fc) + 4r kE ( k ) 
dk 



(fc-fc), 



and then (see the appendix |A"1 for the derivative): 
fcE(fc) w fcE(fc) + |2E(fc) - K(fc)] (fc - fc), 



(34) 



(35) 



which is always bounded since fc < 1. The truncation pro- 
duces an error of the order of (fc — fc) 2 . We can now esti- 
mate the integral of fcE(fc) with respect to z or After 
re-arranging terms, we can write this integral as: 



fcE(fc)dC « Ti(fc)C + T 2 (fc) / fcdC 



where we have defined: 

Ti(fc) = fc [K(fc) - E(fc)] = fc 3 D(fc), 

and 

T 2 (fc) = 2E(fc) - K(fc) = E(fc) - fc 2 D(fc). 



(36) 



(37) 



(38) 



where D(fc) is another common complete elliptic inte- 
gral (see the Appendix [X). These two functions are plotted 
versus fc in Fig. [3] Note that Ti and T 2 logarithmically di- 
verge as fc — > 1. This is not a problem since these functions 
are always invoked with a modulus less than one. Besides, 



7 This is true except when z+ = z— and Z = 0. This case is met 
when the potential is requested in the equatorial plane, at radii 
where the thickness is zero. This is precisely the case of "edges" 
where the mass density p is also expected to vanish. So, we can 
conclude that k < 1 for realistic configurations. 
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2E(k)-K(k) or TJk) 



K'Dfkjoryk) 



0.2 




y point C 



0.4 0.6 
modulus k 



Figure 3. Variation of Ti(fc) and T2(fc) with the modulus k. 
These two functions logarithmically diverge as fc — > 1. 



the integration of fc with respect to £ is easily found analyt- 
ically: 



C v/(a + i?) 2 + C 2 
C 



2v a J? argsh 



(39) 



a + 7? 



By replacing in Eq. (ll8|l the vertical integral by Eq. (|36|l . 
its approximation, we find the "approximate mean Green 
function" : 




(40) 



argsh 



C- 



a + R 



where fc is given by Eq. (|33l) . The error is expected to be of 
the order of f (fc — kydz. The approximation for the po- 
tential is given by Eq.Q where Q is to just be replaced by 
Sapp., or: 



ip(R, Z)^~G E(a)g app . (R, Z; a, z ± ) da. 



(41) 



As this approximation works in the whole physical space 
(no hypothesis has been made upon R and Z) , it can be used 
not only to model the internal structure of self-gravitating 
discs and rings, but also for dynamical studies where the sys- 
tem is basically a source of gravity act i ng on external, mov - 
ing particles (e.g. ISubr fc Karasll2005l ; iTresaco et aHl201lT ). 
Equation (|41[) is also attractive from a numerical point of 
view as it involves a one dimensional integral, with finite 
integrand. 

The long-range properties of the potential are automat- 
ically conserved here as, in deriving this approximation, we 
have made no hypothesis about the relative distance to the 
system. When r — > oo, we have fc ~ 2\faR/r — > 0, and so 
2T 2 (fc) ~ 7T while 4Ti(fc) ~ yrfc 2 (the H- function brings a 
similar, second-order correction; see the Appendix [E)l. The 



axial symmetry 



0.5 - 



point B 




a.R 



Figure 4. A test-torus with circular cross-section, and edges at 
5 and 1 (see Fig. [5] for the mean Green function at points A, B 
and C). 



main contribution in the potential then comes from T2(fc): 



ip(R,Z)^~G I E(a)^/-(— Kk)(-2h)da (42) 



-4ttG- / T,(a)ada 
GM 



8 A NUMERICAL EXAMPLE 

To check the quality of the approximate potential, it is 
not necessary to compare Eq.(|4]) or Eq. (|22[) together with 
Eq. (|41[) . It is sufficient to compare the two Green functions, 
i.e. Q and its approximation <5a PP .. Since the approximation 
made does not involve the H-function (it is the same for 
both expressions), we can simply consider the two members 
of Eq, (|36|) . To perform this comparison, we have considered 
a homogeneous torus with circular cross-section, inner edge 
ai n = | and outer edge a out = 1 as depi cted in Fig. [H A 
simila r test-torus has been considered in iBannikova et al.1 
|201 if ). Also, in order to show how the mean Green function 
behaves (Q also depends on R and Z), we have selected three 
points where the potential could be requested: 

• point A(4,0) inside the torus (at the center), 

• point B(l, 1), outside it, 

• point C (10, 10), relatively "far away" from the system. 

Figure [5] displays the mean Green function Q versus the 
radius a at points A, B and C. As expected (see Section 
01 , the function is the largest when (R, Z) stands inside the 
disc (the case of point A) and peaks at a — R, with a jump 
in the derivative (due to the H-function; see the Appendix 
[C]). At point A, we have fc± w 0.9864, K(fc±) w 3.21 and 
/ kE{k)dz w 0.450 (see i.e. Eq.fl39)), and so 5max. « 8.32 
from Eq. (|20[) . Outside the system (points B and C), Q has 
a much lower amplitude, but still exhibits a minimum. The 
relative difference between Q and its approximation C?a PP . 
is shown in Fig. [6] for points A and B, and in Fig. [7] for 
point C. According to the Taylor expansion (see Sect. [7]), the 
absolute error is given by J (fc — k) 2 dz. This term is easy to 
calculate. Inside the system and in its close neighbourhood, 
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total jump 2it/h(R) 




0.7 O.f 
radius a 

Figure 5. Variation of the mean Green function Q with the radius 
a at points A, B and C for the test-torus considered in Sect. [8] 
(see also Fig. UJ. Inside the system, the radial derivative of the 
5-function undergoes a jump at R = a where the function peaks 

(the case of point A), due to the H-function (see text). 



0.003 



0.0025 



0.002 



<U 

=o 0.0015 



■2 0.001 



0.0005 




Figure 6. Relative difference between the mean Green function 
Q and its approximation C/ a pp. versus a at points A (plain line) 
and B (crosses) for the test-torus shown in Fig. l4l 



6e-06 




0.7 0.8 
radius a 

Figure 7. Same legend as Fig. [5] but for point C. 




-0.4 



0.6 



-0.8 



-1.2 



-1.4 



-1.6 



Figure 8. Potential of the homogeneous torus computed from 
Eq.d22t at the nodes of a regular (R, Z)-grid (64 X 64 nodes; R 
and Z £ [0, 1.5]). The torus has circular cross-section (shown in 
white), a diameter i, and outer edge a ou t = 1. 



exact vs. approximation 




1 1.2 1.4 



Figure 9. Relative error (decimal log. scale) in potential values 
when using the approximation Eq. I|41jl rather than Eq, l|22[l . The 
conditions are the same as for Fig. \E\ 



we find that the error is ~ h 4 /(a + R) 4 . This is a relatively 
small error for geometrically thin discs (most quantities have 
a precision of ~ {h/a) 2 in general). In the present example 
where h reaches 0.25, the absolute error is expected to be 
~ 10 -3 which is compatible with what we observe for points 
A and B. The error in potential values (the area under the 
curves) is therefore of the order of ~ h 4 /(a + R) 4 x Aa 
here, where Aa is the radial extension of the system. This 
is therefore ~ 3 x 10 -4 for point A and for point B, and 
~ 3 x 1CP 6 for point C (see the error map below). 

We have computed the potential map in the vicinity of 
the torus from Eq. Q , using the exact mean Green function 



© ??? RAS, MNRAS 000, [T]-?? 



8 A. Trova, J.-M. Hure and F. Hersant 



Q. The result is shown in Fig. [8] We have also determined 
the approximate potential from Eq. (|41[) . that is, using the 
approximate Green function C/a PP .- The relative error (log. 
scale) between the two potential maps is shown in Fig. [9] 
We conclude that we have, through the approximation, still 
a very precise estimate of J fcE(fc)ri£ inside as well as outside 
the system. This is not that surprising since E(fc) is a regular 
function. 

It is worth noting that theactual test-torus, with 
- ~ 0.33, does not correspond to what is usually called a 
"geometrical thin" system. But as the error is of the order 
0{h 4 /a 4 ), Eq.gU) remains valid for relatively thick systems, 
including thick discs and tori. 



9 INFLUENCE OF VERTICAL 
STRATIFICATION 

As realistic systems are expected to have density gradients in 
the three directions, it is important to check the sensitivity 
of the results obtained so far to the p(z) profile. This can 
easily be done, still assuming axial symmetry, by introducing 
the inhomogeneous analog of Eq.([SJ) defined by: 



gin. 



-p(z)kK(k)dz / / p(z)dz 



(43) 



For geometrically thin disc models, the mass density often 
consists in a power-law of the radius combined with a Gaus- 
sian profile in the direction perpendicular to the midplane, 
which is a direct conseq uence of a loca l ly isothermal gas i n 
hydrostatic equilibrium (jPringld Il98ll : iMiirier et alJl2012h . 
As long as we can decouple the radial and vertical density 
structures, an interesting profile is therefore the following 



/2 / 2z 2 



(44) 



where z 6] — oo, oo[ and po is a constant, possibly a function 
of the radius. We have also considered the parabolic profile, 
namely 



1 - 



(45) 



where z £ [—h, h]. This profile is very similar to the Gaus- 
sian law, matter being however confined into a finite domain. 
For these two inhomogeneous profiles, the surface density is 
2poh, as for the vertically homogeneous case, while central 
values are different (1.5 and « 1.59 and 1, respectively). 
Figure [lOl shows Q ln (a) computed from Eq. (|43[) for the two 
inhomogeneous profiles for the torus considered in the pre- 
vious section, at the interior point A. We see that the mean 
Green functions are very similar in shape, still peaking at 
a — R. The maximum is higher for the parabolic case and 
even higher for the Gaussian case (10% and 11% respectively 
with respect to the canonical case), due to larger and larger 
central valuefl The difference on potential values (the area 



More generally, for a vertical profile of the form: 

, 2 



A+B U 



(46) 



we can find a good approximation for £/ max . by using the right 
expansion for K(fc) (see note [Til followed by integrations terms by 



o 7 



Gaussian profile - 




parabolic profile 



homogeneous case 
(see Fig. 5) 



0.5 



0.6 



0.7 0.S 
radius a 



0.9 



Figure 10. Variation of the mean Green function with the radius 
a at point A for the test-torus considered in Section [8] (see also 
Fig. [JJ for three different profiles p(z) : the homogeneous case 
(see Fig.O, the parabolic profile and the Gaussian profile. 



under the curves) is however not that large, of the order of 
5% (we have w —1.504 and —1.513 respectively, com- 
pared to —1.436 in the homogeneous case). The deviation is 
in excess, as expected: in general, the more concentrated the 
mass distribution, the deeper the potential well. We have 
performed the same comparison at points B and C which 
stand outside the torus. We have noticed that as we move 
away from the distribution, the relative difference between 
profiles gets smaller and smaller. Again, this is not really a 
surprise since at large distance, the potential tends to that 
of a point mass and the details of the distribution are no 
longer perceptible. Figure [TT] shows the relative difference 
between the potential generated by a Gaussian profile and 
the homogeneous case. It turns out that the potential inside 
as well as outside the body is mainly sensitive to the surface 
density, which is therefore the critical parameter. 



10 CONCLUSION 

In this article, we have proposed an alternative definition of 
the potential of axially symmetrical bodies which avoids a 
singular kernel. This is based on the effective integration of 
the genuine kernel along the vertical direction. This natu- 
rally leaves a regularized kernel, called "mean Green func- 
tion" , which peaks at the location of the initial logarithmic 
singularity but remains of finite amplitude. This is yet an- 
other proof that the local contribution of matter located at 
r* « r is very important and may dominate the total contri- 
bution. This new but equivalent form is therefore attractive 
because the corresponding gravitational potential is then ac- 
cessible from two classical integrals: the one over the body's 



terms. We then get: 

5max. 



2Eln 



8R 



2Ah- 



-Bh 



(47) 



h(R) ' 9~ 

which leads to the correct values in the cases considered here, 
namely: ~ 8.32 for the homogeneous profile, ~ 9.02 for the 
parabolic profile and ~ 9.19 for the Gaussian profile. 
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Gaussian vs. homogeneous 




0.2 0.4 0.6 0.1 
R 



Figure 11. Relative error (decimal log. scale) in potential values 
between the Gaussian case and the homogeneous case (see text 
for more details). The conditions are the same as for Fig. [8] 



meridional cross-section, and the other over the body's equa- 
torial radius. The total absence of diverging kernels means 
that there is no need for special meshing or specific numer- 
ical schemes. Also, the properties of Newton's law are fully 
conserved with this approach. 

Rigorously, our results are valid only for mass density 
profiles which are uniform along the z-axis, while there is no 
restriction about the (radial) surface density profile, thick- 
ness/shape and size of the system. This still makes the 
method relevant to various kinds of configurations, rang- 
ing from geome trically thin discs and rings found in dif- 
feren t contexts (ICollin-Souffrin fe Dum ont 199(j; iDubrullel 
1 19921 ; lArnaboldi fe Sparkej 1 1994 ) . to possibly long, cylin- 
drical/filame ntary structures as observed in the interstel- 
lar medium dCurrvl 1200(1 iHennebelld 120031 ) . The fact that 
the asymptotic, long-range behavior of the potential is well 
reproduced means that our formula can also be used to gen- 
erate quickly and simply grids of potential in which the mo- 
tion of test-particles can be studied. Nevertheless, we have 
analyzed a few inhomogeneous situations, in particular the 
Gaussian profile. It turns out that the sensitivity to stratifi- 
cation is weak on potential values. The local surface density 
remains the decisive parameter. 

We failed to convert the mean Green function into 
a line integral — s ee Eq. (|27[) and note [6] — as done in 
lAnsorg et all (|2003h with the genuine kernel fcK(fc). This 
does not seem possible unless new closed-form relationships 
between complete elliptic integrals with homothetical 
moduli are derived. Also, we have expanded at first-order 
the mean Green function around vanishing aspect ratio 
h/a, enabling to perform this conversion. We have then 
obtained an approximation with very small error. Accuracy 
can be improved by considering next terms. 
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APPENDIX A: DEFINITIONS 

The complete elliptic integral of the first, second and third 
kinds are defined by l|Gradshtevn fc Rvzhild[l9o5h : 



K 



(fc) = / - 

Jo J I - fc 2 



rn/2 
E(fc) = / 

Jo 



1 — k 2 sin 2 



and, 



II(m,fc) 



r/2 



(l — m 2 sin 2 <j>) *J\~k 2 sin 2 



(Al) 



(A2) 



(A3) 



respectively, where k is the modulus and m is the charac- 
teristic or parameter. In the present study, we have ^ k ^ 
m ^ 1, k is defined by Eq.©, so that m — k for ( — 0. The 
function D(fe) is defined by: 



fc 2 D(fc) = K(fc) -E(fc). 



(A4) 



The demonstration given in Section[3]is based upon the 
following partial derivatives with respect to the modulus k 
l|Gradshtevn fc RvzhildlT965h : 



%K(k) = - 



E(fc) 



K 



kd k E(k) = E(fe) - K(fc), 
and (|Durandlll953l ): 



9feII(m, k) 



U{m,k) 



E(fc) 



k' 



(A5) 



(A6) 



(A7) 



where k! = y/1 — k 2 is the complementary modulus. 
We also have 



c 



= ± 



V m? — k 2 



2VaR rnk 
where sign + stands for £ > 0, and then 

(drk = -A? (fc 2 - m 2 ) . 
m 2 v ' 



(A8) 



(A9) 




Figure Bl. The function (H(m,mi) normalized to its value at 
x = 0, versus x for different values of m £ [0, 1] . 



APPENDIX B: THE FUNCTION H 

The function H(m, fc) defined in Section [3] can be rewritten 
is terms of the variable x = — £ [0, 1], namely: 

H(m, mx) = mx ^K(mx) — m' 2 II(m, ma:)j , (Bl) 



and so: 



H(m, mx) 2 ,K(mx) — m' H(m,mx) 



H(m,0) 7T 



1 — m' 



where 



H(m,0) = |(l-m'), 



(B2) 



(B3) 



and x' = y/1 — x 2 . The ratio H(m, mx)/H(m, 0) is plot- 
ted versus x and for different parameters m in Fig. IB 11 In 
a first approximation, we have from the figure H(m, fc) w 
x'H(m, 0), and so 

/2 

CH(m,k)tt±Z(a + R)(l-m')—. (B4) 
2 x 



APPENDIX C: JUMP IN THE RADIAL 
DERIVATIVE OF THE FUNCTION H(M, K) 

Whe n a — > R (i.e. m — ¥ 1) and for £ 7^ 0, we have (jPurandl 
Il953h : 



As we have 



7n' 2 n(m, fc) 



d a m' = ± 



-n 2 fc' 
2R 



(a + R) 



2 • 



(CI) 



(C2) 



where sign + stands for a > R, we see that the jump in the 
derivative of n at a — R produces a jump in the derivative 
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of H. This jump is calculated as follows: 

lim d a m' 2 Tl(m, k) = d a lim m' 2 Yl(m, fc) 



(C3) 



7r m 
= ±*- 1 



7r m 1 



2R 



2 k' (a + R) 2 



_ ± 7T 1 1 

2ki2R 



/.From Eg. (|A8[) . we have lim m _n C, = ±2R^. Consequently, 
we get for £ < 0: 



lim <^kd a m' 2 Il(m, k) 

m— >1 



(C4) 



meaning that the jump in the derivative of £H amounts 
to n when a increases and crosses the radius a = R. 
Since the Green kernel contains the difference f+II(m, fc+) — 
£_II(m, fc-), the jump is in fact 2n in total (see for instance 
Fig. [5] point A), except for a point located on the boundary 
(dQ) (i.e. if Z — z±(a)), then it is only n. This is always 
true except for £ = 0. In this case, the jump disappears and 
we have: 



rn'H(m, m) 
which is perfectly continuous. 



E(m) 



(C5) 



APPENDIX D: LONG-RANGE BEHAVIOR FOR 
THE EXACT POTENTIAL 



At large distance from the body, we have 
2VaR 



k < 



0. 



(Dl) 



and so the behavior of t he elliptic integrals is easily d educed. 
In particular, we have (Gr adshtevn fc Rvz hik 1965): 



K(fc) ~ \ 

E(fc) ~ \ 

Tl(m,k) ~ 
We then get: 



(^) 



(D2) 



)] 



-2G 



a)da x j— I kE(k)d( 



-4nG 

GM 
r 



p(a)h(a)ada 
(D3) 



which is the expected result. It can be checked that the con- 
tribution due to the H-function is much smaller and behaves 
like fc 3 . Actually, we find: 



C±H(m,fc±) ~ C±k± — 



fc± 1 ( i k 2 



C±k±j (1 + m 2 



m ) +m 



(D4) 



As k± et m are of the same order, £±H(m, fc±) behaves like 
kl. 



APPENDIX E: LONG-RANGE BEHAVIOR FOR 
THE APPROXIMATE POTENTIAL 



The modulus k± can be approximated as follows: 

k ~ ( x _ aR + Zz ± \ 

r \ r 2 ) ' 

therefore the mean value fc becomes: 



(El) 



2VaR 



2VaR 



1 - 



aR + (z+ + Z-)Z 



aR 



We then get for Ti: 

Tr(fc)(C + -C-) 



(E2) 



(E3) 



Anh- 



laR 



1 - 



aR 



Regarding the second term, we have: 



T 2 (fc) 



h 2 

4 fc 



and 



J kdC, ~ fcC- 

This order is sufficient. We then find 
T 2 (fc) f kdC, 



lit- 1 1 \ 1 — ^fc 2 



(E4) 



(E5) 



(E6) 



and we see that this second term dominates over the first 
one. We then get: 

-2G J p(a)daJ~^ J kE(k)d( 4nG^ J p(a)h(a)ada 

(E7) 



GM 
r 
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